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LONG-TERM  GOALS 

The  long-term  goal  of  our  project  is  to  develop  and  enhance  research  and  educational  capabilities  in 
the  area  of  coastal  engineering  and  science  at  Louisiana  State  University  (LSU)  while  simultaneously 
supporting  the  Navy’s  research  goals  in  the  area  of  Coastal  Geosciences.  The  focus  of  the  present  work 
is  to  develop  a  new  modeling  framework  for  simulations  of  coastal  processes  in  deltaic  enviromnents 
using  advanced  numerical  methods  and  high  performance  computing  technology.  In  particular,  the 
utilization  of  adaptive  numerical  methods  such  as  the  spectral  element  method  on  modern  computer 
platforms  with  thousands  of  multi-core  processors  will  enable  coastal  modelers  to  simulate  complex 
physical  processes  with  improved  accuracy  and  efficiency. 


OBJECTIVES 

The  specific  objectives  of  this  project  are  to: 
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•  Develop  the  capability  of  modeling  coastal  circulation  and  nearshore  surface  waves  in  deltaic 
sedimentary  and  hydrodynamic  environments  in  an  integrated  modeling  framework  by 
extending  the  Boussinesq  theory  for  nearshore  hydrodynamics  to  muddy  coasts  and  non¬ 
hydrostatic  three-dimensional  (3D)  flow  regimes  with  stratifications. 

•  Apply  advanced  high-order  numerical  methods  to  coastal  hydrodynamic  modeling. 

•  Simulate  large-scale,  long-tenn  problems  in  the  deltaic  environment  by  integrating  the 
application-oriented  modeling  system  with  massive-processor  computing  facilities  and 
technologies  available  at  LSU. 

•  Complement  the  Office  of  Naval  Research  recent  research  initiatives  on  Tidal  Flats  and  Wave- 
Mud  Interactions  by  integrating  the  new  modeling  system  with  the  field  data  collected  in  those 
programs. 


APPROACH 

The  research  project  consists  of  theoretical  formulation  and  analysis,  the  development  and  verification 
of  an  advanced  modeling  system  using  the  spectral/hp  element  methods  and  high-perfonnance 
computing  technologies,  and  the  utilization  of  the  new  model  as  a  research  tool  to  advance  knowledge 
and  understanding  of  coastal  circulation  and  nearshore  waves  in  deltaic  sedimentary  and  hydrodynamic 
environments.  Relevant  methodologies  in  three  disciplines:  civil  engineering,  physical  oceanography 
and  computational  science,  are  being  utilized.  Interdisciplinary  interactions  are  taking  place  among  the 
investigators  in  different  fields  and  through  recruiting  graduate  students  from  the  three  disciplines.  All 
project  members  meet  together  weekly,  with  meetings  alternating  between  focusing  on  coastal  science 
and  computational  science.  Code  development  and  document  preparation  is  enabled  through  a  project¬ 
wide  source  code  versioning  system. 

A  new  approach  is  taken  to  meet  the  objectives  of  this  project:  1)  Use  the  Boussinesq  theory  (e.g., 
Chen  2006)  to  improve  the  efficiency  of  non-hydrostatic  3D  Navier-Stokes  equation  solvers  as  well  as 
to  extend  the  applicability  of  the  modeling  system  to  deltaic  environments,  and  2)  utilize  spectral/hp 
element  methods  with  unstructured  grids  to  solve  the  partial  differential  equations  (PDE)  under 
realistic  deltaic  conditions  on  high-performance,  massive-processor  computers  available  at  LSU. 

Boussinesq-type  models  for  breaking-driven  currents  in  a  large  nearshore  area  are  computationally 
demanding.  Taking  seabed  conditions  into  account  by  the  Boussinesq  model  will  further  increase  the 
computational  effort  by  a  factor  of  two.  It  is  therefore  desired  to  speed  up  Boussinesq  models  for 
practical  applications,  in  particular  for  morphological  and  ecological  simulations.  In  addition,  a  direct 
solution  to  the  Euler  equations  of  motion  for  coastal  waves  and  currents  is  also  desirable  for  fully 
dispersive  waves  and  currents  with  strong  vertical  shears.  The  solution  to  the  growing  demand  of 
computing  power  by  coastal  models  is  the  use  of  high-perfonnance  computing  (HPC)  technologies. 


WORK  COMPLETED 

Wave  attenuation  and  mass  transport  of  a  water-mud  system  due  to  a  solitary  wave  on  the  free  surface 
has  been  modeled  by  using  the  Chebyshev-Chebyshev  collocation  spectral  method  for  spatial 
discretization  and  a  fourth-order  multi-stage  scheme  for  time  integration.  The  governing  equations 
were  formulated  in  Lagrangian  coordinates  and  perturbation  equations  for  shallow  water  waves  were 
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derived.  An  iteration-by-subdomain  technique  was  introduced  to  tackle  the  interface  in  the  two-layer 
system.  The  numerical  model  was  tested  against  available  analytical  solutions  and  good  agreement  has 
been  found.  Although  the  model  is  focused  on  solitary  waves  and  Newtonian  fluid-mud,  the 
methodology  can  be  extended  to  oscillatory,  nonlinear  water  waves  over  a  non-Newtonian  mud 
bottom. 

The  Boussinesq  approach  has  been  applied  to  a  three-dimensional  Euler  solver  for  highly  dispersive 
waves  and  stratified  flows.  Free  from  the  irrotational  flow  assumption,  the  three-layer  model  can 
correctly  resolve  the  free-surface  displacement  and  the  velocity  profile  of  short  waves  interacting  with 
shear  currents,  which  cannot  be  easily  obtained  by  the  existing  Boussinesq-type  equations. 

Considerable  effort  has  been  devoted  to  the  implementation  of  the  discontinuous  Galerkin  scheme 
using  the  open-source  spectral//?/;  library  Nektar++  (www.nektar.info)  and  its  predecessor,  Nektar. 
Through  Nektar++/Nektar,  the  fundamental  routines  associated  with  a  high-order  finite-element 
method  are  easily  accessible.  We  have  solved  the  nonlinear  shallow  water  equation,  the  Boussinesq- 
type  equations  for  waves  and  currents,  and  the  Euler  equations  of  motion  on  the  vertical  plane  based  on 
spectral  element/hp  discontinuous  Galerkin  methods.  Testing  of  the  codes  is  being  conducted. 

Our  approach  in  parallelizing  the  shallow  water  equation  (SWE)  solver  was  to  integrate  our  code  into 
the  Cactus  computational  framework  (Goodale  et  al.  2003)  and  provide  parallelism  through  a  Cactus 
“Driver”  module  (or  “thorn”).  Cactus  (www.cactuscode.org)  is  an  open  source  problem  solving 
environment  designed  for  scientists  and  engineers,  its  modular  structure  easily  enables  parallel 
computation  across  different  architectures  and  collaborative  code  development  between  different 
groups.  Cactus  originated  in  the  academic  research  community,  where  it  was  developed  and  used  over 
many  years  by  a  large  international  collaboration  of  physicists  and  computational  scientists. 

In  order  to  integrate  the  serial  SWE  solver  into  Cactus,  several  thorns  have  been  developed.  Most 
importantly,  we  have  designed  a  “Nektar++”  thorn  that  initializes  and  populates  the  data  structures  of 
Nektar++.  We  also  provide  a  “SWE”  thorn  that  contains  the  actual  SWE  solver  based  on  routines 
defined  in  the  SWE  solver  library.  Both  spectral/ph  element  SWE  and  fourth-order  finite  difference 
SWE  solvers  have  been  integrated  with  Cactus  to  take  advantage  of  its  parallelism. 


RESULTS 

The  major  results  obtained  in  2009  are:  1)  the  development  of  a  high-order  numerical  model  for  wave 
attenuation  and  mass  transport  of  a  wave-mud  system  due  to  a  solitary  wave  using  the  Chebyshev- 
Chebyshev  collocation  spectral  method  for  spatial  discretization  and  a  fourth-order  multi-stage  scheme 
for  time  integration,  2)  the  improvement  of  a  three-dimensional  (3D)  Euler  solver  for  fully  dispersive 
water  waves  by  extending  the  Boussinesq  approach  to  the  treatment  of  the  non-hydrostatic  pressure  on 
the  top  layer  of  the  3D  model,  3)  the  continued  implementation  and  testing  of  the  Boussinesq  model 
(FUNWAVE)  for  waves  and  currents  on  a  solid  bed  into  the  Cactus  Framework  using  high-order 
finite-difference  schemes  to  serve  as  a  verification  tool,  4)  the  implementation  of  spectral/hp  element 
method  for  incompressible  flow  with  and  without  a  free  surface,  and  5)  the  public  release  of  the 
Shallow  Water  Solver  as  part  of  the  Nektar  ++  software  package  for  computational  fluid  dynamics 
based  on  the  spectral/hp  element  method  (www.nektar.info). 
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A  Spectral  Collocation  Model  for  Solitary  Wave  Attenuation  and  Mass  Transport  over  Viscous  Mud 

The  numerical  model  solves  the  first-  and  second-order  equations  derived  from  the  Navier-Stokes 
equations  for  viscous  fluid  motion  in  Lagrangian  coordinates.  Shallow  water,  week  displacement,  and 
a  horizontal  bed  have  been  assumed  in  the  perturbation  analysis  for  a  two-layer  Newtonian  fluid-mud 
system  on  the  vertical  plane.  The  numerical  model  based  on  the  high-order  numerical  methods  has 
been  tested  against  available  analytical  solutions  and  good  agreement  has  been  found.  Due  to  the  high 
viscosity  of  the  bottom  fluid-mud,  the  modeled  horizontal  particle  velocities  and  particle  displacements 
decrease  with  the  distance  as  the  wave  propagates,  which  justifies  the  use  of  non-periodic  spectral 
schemes.  The  computed  wave  attenuation  rate  is  in  good  agreement  with  the  analytical  solution  of 
Jiang  and  Zhao  (1989)  when  the  thickness  of  the  mud  layer  is  relatively  thick.  Numerical  simulations 
however  suggest  that  the  accuracy  of  the  existing  boundary  layer  theory  for  wave-mud  interaction  is 
limited  when  the  mud-layer  thickness  becomes  smaller  than  the  wave  height  because  the  assumed 
irrotational  core  may  not  exist.  Although  the  model  has  been  focused  on  solitary  waves  and  Newtonian 
fluid-mud,  the  modeling  framework  developed  can  be  extended  to  oscillatory,  nonlinear  waves  over  a 
non-Newtonian  mud  bottom  modeled  as  a  power-law  fluid.  The  result  was  published  in  the  Journal  of 
Engineering  Mechanics  (Huang  and  Chen  2009). 

Improvement  of  an  Euler  Solver  for  Fully  Dispersive  Water  Waves  Using  the  Boussinesq  Approach 

The  basic  idea  to  improve  the  accuracy  of  a  three-dimensional  Euler  solver  for  fully  dispersive  waves 
is  to  obtain  an  analytical-based  fonn  of  non-hydrostatic  pressure  distribution  at  the  top  layer  by  using 
the  principle  of  Boussinesq  theory:  1)  The  vertical  profiles  of  horizontal  velocity  components  are 
obtained  by  using  the  Taylor  series  expansion  at  the  reference  location;  2)  the  vertical  velocity  profile 
is  obtained  by  substituting  the  horizontal  velocity  components  into  the  continuity  equation  and 
integrating  the  resulting  equation  from  the  reference  to  an  arbitrary  location  z  ;  and  3)  the  analytical 
form  of  pressure  profile  is  obtained  by  integrating  the  vertical  momentum  equation  from  an  arbitrary  z 
to  the  free  surface  and  applying  the  Leibniz’s  rule  with  the  use  of  free-surface  kinematic  boundary 
condition  and  free-surface  pressure  condition.  Numerical  experiments  have  shown  that  this 
Boussinesq-like  treatment  of  the  top  layer  enables  the  use  of  as  few  as  four  vertical  layers  to 
accurately  model  nonlinear  short  waves  with  dimensionless  wave  number,  kh=16,  which  can  not  be 
easily  simulated  by  existing  Boussinesq-type  models.  Different  from  conventional  Boussinesq  models, 
the  integration  in  the  vertical  direction  is  only  limited  to  the  top  layer  and  no  assumption  of  zero 
horizontal  vorticity  components  was  made.  The  result  will  appear  in  the  Journal  of  Waterway,  Port, 
Coastal  and  Ocean  Engineering  (Wu  et  al.  2010). 

A  Highly  Parallelized  Shallow  Water  Wave  Solver  Based  on  Combined  OpenMP  and  MPI 

As  part  of  the  continued  development  of  the  parallel  Boussinesq  model  for  waves  and  currents  within 
the  Cactus  Framework  using  high-order  finite-difference  schemes,  a  combined  OpenMP  and  MPI  have 
been  used  to  speed  up  the  code.  We  used  4th  order  central  differencing  for  spatial  discretization  and  4th 
order  Runge-Kutta  method  for  time  integration.  A  ShallowWater  thorn  (or  module)  implemented  in 
Cactus  uses  Carpet  as  the  parallel  driver.  Figure  1  shows  the  result  of  weak  scaling  test  of  the  shallow 
water  wave  solver  for  wave  sloshing  in  a  rectangular  tank.  The  number  of  processors  ranges  from  4 
cores  (1  node)  to  256  cores.  The  mesh  size  was  1004x1004  (including  the  grid  points  used  for 
imposing  the  wall  boundary  conditions)  in  the  x-  and  y-direction  for  the  4  core  job  and  increased 
proportionally  as  the  number  of  cores  increases.  The  computational  work  per  core  is  maintained 
constant  and  the  total  size  of  the  problem  increases  with  increase  in  number  of  processors  (weak 
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scaling).  For  an  ideal  weak  scaling,  the  time  to  the  solution  should  not  increase  significantly  when  the 
problem  size  is  increased  in  proportion  to  the  number  of  processors.  All  these  tests  were  performed  on 
the  “Painter”  supercomputer  (www.loni.org)  that  has  680  nodes  with  each  node  having  2  quad-core 
processor  (i.e.  8  cores)  with  8  Gb  RAM.  It  is  seen  in  Figure  1  that  good  scaling  has  been  achieved  up  to 
256  cores  which  was  the  limit  of  these  initial  experiments.  Since  the  Cactus  parallel  drivers  have 
already  been  shown  to  scale  to  over  130,000  cores  with  single  grid  finite  difference  problems  we 
expect  this  scaling  to  continue.  The  variance  in  the  scaling  seen  in  the  graph  is  likely  due  to  the 
changing  topology  of  the  decomposition  of  the  computational  grid  across  the  cores  which  causes 
variations  in  the  local  grid  sizes  on  each  core.  The  different  lines  on  the  graph  relate  to  the  different 
numbers  of  OpenMP  threads  used. 


#  cores 

Figure  1  Weak  scaling  tests  of  the  shallow  water  wave  solver  in  Cactus. 


Continuous  and  Discontinuous  Galerkin  Spectral/hp  Element  Methods  for  Incompressible  Flows 

For  complex  3D  flow  problems  with  strong  vertical  variations,  a  direct  solution  to  Euler  equation  or 
Navier-Stokes  (N-S)  equation  is  desirable.  Basically,  there  are  two  approaches  for  solving  the  pressure 
field  applied  to  incompressible  flow  problems  in  an  Euler  solver.  The  first  method  is  to  solve  the 
Helmholtz  or  Possion  equation  to  obtain  the  pressure  field  (e.g.,  Karniadakis  and  Sherwin  2005)  and 
the  second  one  is  the  artificial  compressibility  (AC)  method  first  proposed  by  Chorin  (1967).  The 
attributes  of  AC  methods  are  its  high  efficiency  (less  memory  requirement  compared  to  directly 
solving  Helmholtz  equation),  parallel  computing,  and  its  robustness  for  complicated  flow  problems. 
The  implementation  of  AC  method  is  relative  simple.  Comparisons  of  both  approaches  in  the  context 
of  discontinuous  Galerkin  spectral/hp  element  methods  for  water  waves  are  being  carried  out. 

In  order  to  simulate  nonlinear  water  waves,  the  arbitrary  Lagrangian-Eulerian  (ALE)  method  is  utilized 
to  handle  the  moving  free  surface  in  the  current  spectral  element  code.  ALE  is  a  well-documented 
algorithm  for  moving  domain  and  deformation  problems  (Hirt  et  al.  1974).  We  are  benchmarking  the 
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code  to  determine  the  performance  of  the  ALE  for  highly  nonlinear  water  wave  problems  solved  by 
spectral/hp  element  methods. 

The  spectral/hp  element  solver  for  the  N-S  equation  is  being  tested  and  the  preliminary  results  are 
shown  in  Figures  2-6.  The  first  test  case  is  a  backward-facing  flow  (channel  flow)  in  which  the  step  is 
a  half  height  of  the  channel  width  (Figure  2).  As  observed  in  many  experiments,  the  steady  discharge 
in  the  upstream  produces  a  large  re-circulation  zone  attached  to  the  backward-facing  step,  which  is  a 
useful  test  case  for  numerical  models.  Figures  3-4  illustrate  the  computed  flow  separation  and  the 
resultant  velocity  and  pressure  fields. 


Figure  2  Schematic  of  a  backward  facing  flow  setup  for  benchmarking  the  spectral/hp  element  model. 


Figure  3  Modeled  velocity  vectors  of  the  backward-facing  flow  at  Re- 1000  using  quadrilateral 
elements  at  a  transition  state  (t=4  s).  The  inlet  velocity  is  u=  2.0  m/s  and  v=0.0  m/s,  and  the 

downstream  boundary  is  non-reflective. 
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Figure  4  Pressure  contours  for  a  backward  facing  flow  at  Re= 1000  and  t=4s  as  in  Figure  3. 


The  second  test  case  is  an  incompressible  flow  passing  a  horizontal  circular  cylinder.  Unstructured 
triangular  elements  are  applied  in  the  computational  domain  with  curved  edges/surfaces.  Numerical 
results  from  the  spectral/hp  element  code  show  that  the  vertex  shedding  starts  at  Re=116,  which  is  in 
very  good  agreement  with  other  finite  element  simulations  with  much  finer  spatial  resolution  (e.g., 
Kundu  and  Cohen  2002).  More  quantitative  comparisons  between  the  results  from  the  new  solvers  and 
available  analytic  solution  and  experimental  data  are  underway. 


Figure  5  Computed  Vortex  shedding  by  the  spectral/hp  element  method  with  unstructured  triangular 
elements  at  Re=150.  The  shedding  starts  at  Re=116,  which  is  very  close  to  the  result  of  a  body-fitted 
finite  element  simulation  with  finer  resolution  (Kundu  and  Cohen,  2002). 
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Figure  6  Computed  pressure  contours  of  the  flow  passing  a  horizontal  cylinder  at  Re=150. 


Public  Release  of  the  Discontinuous  Galerkin  Shallow  Water  Solver  in  Nektar  ++ 

The  close  collaboration  between  the  LSU’s  team  and  the  Nektar  ++  developers  (Prof.  Spenser  Sherwin 
and  his  team)  at  the  Imperial  College  London,  U.K.  has  resulted  in  the  public  release  of  the  Spectral/hp 
element  discontinuous  Galerkin  (DG)  shallow  water  solver  (SWS)  developed  though  this  project  as 
part  of  the  latest  version  of  the  Naktar  ++  software  package  for  computational  fluid  dynamics 
(www.nektar.info).  The  DG  SWS  class  contains  functions  for  the  evaluation  of  the  flux  vectors, 
numerical  fluxes,  equation  dependent  boundary  conditions,  various  source  terms,  etc.  This  class 
provides  a  shallow  water  equation  solver  library  in  Nektar++.  Details  on  the  numerical  implementation 
of  the  DG  SWS  have  been  documented  in  Eskilsson  et  al.  (2009). 


IMPACT/APPLICATIONS 

We  have  successfully  implemented  a  spectral//?/)  discontinuous  Galerkin  method  for  solving  the  SWE 
and  interfaced  it  with  the  Cactus  computational  framework.  The  DG  shallow  water  solver  has  become 
a  part  of  the  Nektar++  package  for  CFD  available  to  the  public.  We  have  also  developed  a  spectral 
model  to  study  wave  and  mud  interaction,  and  have  improved  the  accuracy  of  a  3D  Euler  solver  for 
highly  dispersive  waves.  We  have  obtained  encouraging  results  in  terms  of  parallel  performance  and 
the  coupled  software’s  scaling  ability.  The  use  of  Cactus  provides  a  path  for  extensibility,  integrating 
with  cutting-edge  computational  hardware  and  cyberinfrastructure,  and  building  a  comprehensive 
toolkit  for  coastal  applications.  The  results  were  documented  in  the  recent  publications  (Huang  and 
Chen  2009,  Wu  et  al.  2010,  and  Eskilsson  et  al.  2009). 
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The  research  is  expected  to  improve  the  Navy’s  capability  of  modeling  nearshore  surface  waves  and 
coastal  processes  in  heterogeneous  sedimentary  environments.  The  modeling  framework  integrated 
with  the  CFD  Toolkits  developed  at  LSU  will  allow  us  to  improve  the  accuracy  of  hydrodynamic 
models  and  couple  them  with  sediment  transport  models  for  coastal  morphodynamic  studies. 

In  addition  to  supporting  the  Navy’s  research  goals,  the  proposed  project  is  contributing  to  the 
Louisiana  State  University’s  mission  on  research  and  graduate  education.  The  north  Gulf  Coast,  where 
the  Naval  Research  Laboratory  and  other  naval  facilities  are  located,  is  in  need  of  research  and 
education  in  coastal  engineering.  The  training  of  post-doctoral  fellows  and  graduate  students  will 
enhance  the  graduate  grogram  in  coastal  engineering  at  LSU  to  meet  the  need  for  graduate  education 
on  the  Gulf  Coast  in  support  of  national  defense. 

TRANSITIONS 

The  spectral/hp  element  discontinuous  Galerkin  shallow  water  solver  developed  through  this  project 
became  a  part  of  the  open  source  Nektar++  software  package  for  computational  fluid  dynamics 
available  to  the  public  (www.nektar.info). 

RELATED  PROJECTS 

Our  project  is  leveraging  and  coordinating  with  activities  in  several  other  ongoing  activities: 

XiRel:  This  NSF  funded  project  is  optimizing  and  extending  an  Adaptive  Mesh  Refinement  layer  for 
the  Cactus  framework,  which  will  be  used  for  our  structured  grid  codes. 
(http://www.cactuscode.org/Development/xirel) 

ALPACA:  This  NSF  funded  project  is  developing  debugging  and  profiling  tools  for  the  Cactus 
framework  which  will  support  the  Coastal  Modeling  Framework  developed  in  this  project. 
(http://www.cactuscode.org/Development/alpaca) 

CyberTools:  This  NSF/BOR  funded  project  is  developing  a  cyberinfrastructure  across  the  100  TFlop 
machines  of  the  Louisiana  Optical  Network  Initiative.  Our  project  is  providing  one  of  the  application 
drivers  for  this  infrastructure,  (http://cybertools.loni.org) 

Blue  Waters:  This  large  NSF  funded  project  will  deploy  a  petascale  computing  facility  at  NCSA  in 
2011.  Our  Cactus-based  coastal  application  will  provide  a  driver  problem  for  the  development  of 
software  for  this  resource. 

CFD  IGERT:  An  NSF  graduate  training  and  education  program  at  LSU  in  training  students  in 
computational  fluid  dynamics  and  high  performance  computing.  Several  research  projects  are  building 
on  the  CFD  Toolkit  which  is  contributing  to  our  project. 

SCOOP:  Where  appropriate,  our  models  will  be  integrated  into  the  community  infrastructure  of  the 
NOAA/ONR  funded  SURA  Coastal  Ocean  and  Observing  Program.  SCOOP  maintains  a  coastal 
archive  at  LSU  with  realtime  forcing  and  simulation  data  for  stonn  events. 
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NSF-CAREER:  The  five-year  research  project  is  focused  on  simulations  of  nonlinear  coastal  waves 
and  air-sea  momentum  fluxes,  which  complements  the  present  research  project. 
(http://www.nsf.gov/eng/cbet/nuggets/1443/1443_chen.htm) 
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